function [z mix] = gmmprod(mix1, mix2)

if nargin < 2
    mix2 = mix1;
end

mix = gmm(mix1.nin, mix1.ncentres*mix2.ncentres, 'full');

count = 0;
for i = 1:mix1.ncentres
    for j = 1:mix2.ncentres
        count = count+1;
        
        [zz mix.centres(count, :) mix.covars(:, :, count)] = ...
            gaussprod(mix1.centres(i, :), mix1.covars(:, :, i), mix2.centres(j, :), mix2.covars(:, :, j));
        mix.priors(count) = mix1.priors(i)*mix2.priors(j)*zz;
    end
end

z = sum(mix.priors);
mix.priors = mix.priors/z;
